Ubiquitous lognormal distribution of neuron densities in mammalian cerebral cortex

Abstract Numbers of neurons and their spatial variation are fundamental organizational features of the brain. Despite the large corpus of cytoarchitectonic data available in the literature, the statistical distributions of neuron densities within and across brain areas remain largely uncharacterized. Here, we show that neuron densities are compatible with a lognormal distribution across cortical areas in several mammalian species, and find that this also holds true within cortical areas. A minimal model of noisy cell division, in combination with distributed proliferation times, can account for the coexistence of lognormal distributions within and across cortical areas. Our findings uncover a new organizational principle of cortical cytoarchitecture: the ubiquitous lognormal distribution of neuron densities, which adds to a long list of lognormal variables in the brain.


Introduction
Neurons are not uniformly distributed across the cerebral cortex; their density varies strongly across areas and layers (Brodmann 1909;Economo et al. 2008). The neuron density directly affects short-range as well as long-range neuronal connectivity (Braitenberg and Schüz 1991;Ercsey-Ravasz et al. 2013). Elucidating the distribution of neuron densities across the brain therefore provides insight into its connectivity structure and, ultimately, cognitive function. Additionally, statistical distributions are essential for the construction of computational network models, which rely on predictive relationships and organizational principles where the experimental data are missing (Hilgetag et al. 2019;Albada et al. 2022). Previous quantitative studies have provided reliable estimates for cell densities across the cerebral cortex of rodents (Herculano-Houzel et al. 2013;Charvet et al. 2015;Erö et al. 2018), non-human primates (Collins et al. 2010;Charvet et al. 2015;Collins et al. 2016;Turner et al. 2016;Atapour et al. 2019;Beul and Hilgetag 2019), large carnivores (Jardim-Messeder et al. 2017), and humans (Economo et al. 2008;Bartheld et al. 2016). However, to the best of our knowledge, the univariate distribution of neuron densities across and within cortical areas has not yet been statistically characterized. Instead, most studies focus on qualitative and quantitative comparisons across species, areas, or cortical layers. Capturing the entire distribution is necessary because long-tailed, highly skewed distributions are prevalent in the brain (Buzsáki and Mizuseki 2014) and invalidate the intuition-guided by the central limit theorem-that the vast majority of values are in a small region of a few standard deviations around the mean.
Here, we characterize the distribution of neuron densities ρ across mammalian cerebral cortex. Based on the sample histograms ( Fig. 1) we hypothesize that ρ follows a lognormal distribution, similar to many other neuroanatomical and physiological variables (Buzsáki and Mizuseki 2014) such as synaptic strengths (Robinson et al. 2021), synapse sizes (Loewenstein et al. 2011;Santuy et al. 2018;Dorkenwald et al. 2022), axonal widths (Wang et al. 2008;Liewald et al. 2014), and corticocortical connection densities (Markov et al. 2014;Gȃmȃnuţ et al. 2018). We used neuron density data from mouse (Mus musculus), marmoset (Callithrix jacchus), macaque (Macaca mulatta), human (Homo sapiens), galago (Otolemur garnettii), owl monkey (Aotus nancymaae), and baboon (Papio cynocephalus anubis) to test this hypothesis (see Section Cell density data for a detailed description of the data). The marmoset, galago, owl monkey, baboon, and macaque 2 data sets are based on a single subject; the mouse, macaque 1 , and human data sets are based on a combination of data across several subjects. The statistical tests conclude that the hypothesis cannot be rejected in the majority of cases, suggesting that the underlying distribution is compatible with a lognormal distribution if the samples are based either on cytoarchitectonically defined areas or on uniformly sampled regions. Beyond the distribution across cortical areas, we show that neuron densities within most areas of marmoset cortex are also compatible with a lognormal distribution. To complement the statistical tests, we perform a model comparison with several other distributions and find that none outperform the lognormal distribution as a model of the data within and across areas. Finally, we show that the lognormal distribution within cortical areas can emerge during neurogenesis from a simple cell division model with variability. The model can furthermore account quantitatively for the lognormal density distribution across areas based on an inferred distribution of proliferation times of the Fig. 1. Neuron and cell densities ρ follow a lognormal distribution across cortical areas for multiple species. Different data sets for the same species are denoted with subscript indices (see Section Cell density data). (A) Histogram of ρ (bars) and probability density function of a fitted lognormal distribution (line). The number of samples N is either the number of sampled cytoarchitectonic areas (mouse, marmoset, macaque 1 , human, galago 1 , and owl monkey) or the number of sampling frames (baboon, macaque 1 , and galago 2 ). For marmoset, galago, owl monkey, baboon, and macaque 2 , the data are based on a single subject; for mouse, macaque 1 , and human it is based on a combination of data across several subjects. (B) Z-scored ln(ρ) histogram (bars), standard normal distribution (line), and result of the SW normality test. (C) Probability plot of z-scored ln(ρ). Discarded outliers marked with a red cross.
areas. Additional between-area variability in the proliferation rates is also compatible with the model but not necessary to obtain a quantitative agreement with the observed distribution. Thus, our model shows how the lognormal distribution of neurons could emerge both within and across the cortical areas.
In the cases of mouse, marmoset, macaque 1 , human, galago 1 , and owl monkey the data were sampled from standard cytoarchitectonic parcellations; abbreviated names for all areas are listed in Table S1. Note that we use subscript indices to distinguish between different data sets on the same model animal, e.g. macaque 1 and macaque 2 .
Neuron density estimates for the mouse were published by Erö et al. (2018), and were measured from previously published Nisslbody-stained slices (Lein et al. 2007), where genetic markers were used to distinguish between cell types. The data were provided in the Allen Brain Atlas parcellation (Lein et al. 2007;Dong 2008). To estimate the quantitative neuron and cell densities Erö et al. (2018) used "a variety of whole brain image datasets [from the Allen Mouse Brain Atlas]" as well as "some values reported from anatomical experiments in the literature", such as from (Herculano-Houzel et al. 2013).
Neuron density estimates for the marmoset cortex were published by Atapour et al. (2019), and were measured from NeuNstained slices. The data were provided in the Paxinos parcellation (Paxinos et al. 2011) and all quantitative values are derived from the same brain of a single subject. Neuron densities within each counting frame used in the original publication (Atapour et al. 2019, their Fig. S1) were obtained via personal communication with Nafiseh Atapour, Piotr Majka, and Marcello G. Rosa.
The neuron density estimates in the first macaque data set, macaque 1 , were previously published in visual form by Beul and Hilgetag (2019), and were obtained from both Nissl-body-and NeuN-stained brain slices. Most of the numerical values were reported by Dombrowski et al. (2001) and Hilgetag et al. (2016), and the remaining values were provided by Sarah F. Beul and Claus C. Hilgetag via personal communication. Counts based on Nissl-body staining were scaled according to a linear relationship with the counts from NeuN staining obtained from selected areas where both types of data were available (Beul and Hilgetag 2019). The data follow the M132 parcellation (Markov et al. 2014) and constitute the average across subjects.
Cell density estimates for the human cortex were previously published by Economo et al. (2008), and were measured from Nissl-body-stained brain slices. The human data therefore most likely ref lect combined neuron and glia densities. The data were provided in the von Economo parcellation (Economo et al. 2008); all quantitative values are based on "the mean of the numbers gathered from the various brains" (Economo et al. 2008, p. 201]. Cell and neuron density estimates for galago 1&2 , owl monkey, baboon, and macaque 2 were previously published by Collins et al. (2010), and were measured using the isotropic fractionator method. The data are sampled from common parcellation schemes in galago 1 and owl monkey, approximately equal-size samples in the baboon, and non-uniform samples in macaque 2 and galago 2 . We refer to uniform sampling when the cortex was subdivided into small regions of approximately equal size and shape presumably without large cytoarchitectonic variations (the case for the baboon), and to non-uniform sampling when the samples often cross cytoarchitectonic boundaries (the case for macaque 2 and galago 2 ). All the quantitative values were derived from one hemisphere of each subject separately. Each data set thus describes only a single subject, without combining data across individuals of the same species.
Note that all the samples representing either full cortical areas or crossing area boundaries include the entirety of the gray matter, spanning all layers of cortex. The only layer-resolved data set is the marmoset; here the mean across all layers was reported as the density.

Lognormality testing
To test for lognormality, we take the natural logarithm, ln(ρ), which converts lognormally distributed samples into normally distributed samples. We then test for normality of ln(ρ) using the Shapiro-Wilk (SW) test, the most powerful among a number of commonly applied normality tests (Razali and Yap 2011). Large outliers (|z-scored ln(ρ)| ≥ 3) were excluded from the normality test; the excluded outliers are indeed cytoarchitectonically distinct areas (discussed below).
Note that any hypothesis test, including the SW test, cannot show that the distribution is lognormal. If the p-value is larger than a certain threshold one cannot reject the null hypothesis that the distribution of ρ is lognormal, i.e. the data are compatible with a lognormal distribution. Thus, we perform further tests, such as statistical model comparison, and comparing ln(ρ) across different animals with a Kolmogorov-Smirnov test.

Statistical model comparison
In order to assess which model is most compatible with the data, we compared the relative likelihood of different distributions against each other. We included an extensive list of distributions with support on R + , estimated the distributions' parameters using maximum likelihood, and calculated the Akaike Information Criterion (AIC) where k is the number of estimated parameters of the model and L is the estimated maximum likelihood. We further compare the models using the relative likelihood (L r ) where AIC min is the minimum AIC across all models and AIC i is the AIC for the ith model. The relative likelihood indicates the probability that, from among the tested models, the ith model most strongly limits the information loss (Burnham and Anderson 2004). We take a significance threshold of α = 0.05 on the relative likelihood to determine whether a model is significantly worse than the best possible model.

Within areas
Neurons are generated through symmetric or asymmetric cell division of neural progenitor cells (Cadwell et al. 2019). Grouping all types of progenitor cells into a single population P, all neurons into a population N, and excluding other post-mitotic cell types, their population size can be modeled by the coupled system of differential equations (Picco et al. 2018) where λ P denotes the (potentially time-dependent) rate of progenitor generation and λ N the (potentially time-dependent) rate of neuron generation. Following the radial unit hypothesis (Rakic 2009), we consider a small number of such radial units (small compared to the total size of the area) and determine the density of progenitors ρ P and neurons ρ N in these radial units by dividing through the volume V of the considered radial units in the fully developed cortex. Importantly, this reference volume is the same for every area. Since Equation (3) is linear, dividing by V leads to Note that this does not necessarily exclude tangential migration as long as the net inf luence of the tangential migration is zero.
Progenitor cell proliferation First, we consider the proliferation of the progenitor cells, which we assume to be governed by a noisy rate. Modeling the noise by a zero-mean, unit-strength Gaussian white noise process ξ , we obtain a stochastic differential equation (SDE) for the progenitor cell density where σ controls the (potentially time-dependent) intensity of the noise. Using the Stratonovich interpretation-assuming that the noise process has a small but finite correlation time before taking the white-noise limit (Van Kampen 2007)-the SDE transforms by the same rules as an ordinary differential equation (Van Kampen 2007) such that we can rewrite the SDE (5) Since ξ(t) is Gaussian and Equation (6) is linear, ln ρ P (t) is Gaussian and hence ρ P (t) is lognormally distributed at all times t. The parameters of this lognormal distribution are the mean of the logarithmic progenitor cell density μ P (t) = ln ρ P (t) and the variance of the logarithmic progenitor cell density σ P (t) 2 = (ln ρ P (t)) 2 (here x ≡ x − x ). Using Equation (6), ξ(s) = 0, and ξ(s)ξ(s ) = δ(s − s ), we obtain (cf. for instance Braumann (2007)) Thus, the progenitor cell densities are lognormally distributed at all times with parameters μ P (t) and σ P (t) 2 . The corresponding first 2 moments of the progenitor cell density are where we used the characteristic functional of the Gaussian white noise Neurogenesis Next, we consider neurogenesis. We assume that the noise affects primarily the rate of progenitor cell proliferation. The solution to Equation (4) for a given ρ P (t) and initial condition Since ρ N (t) is the integral of the (temporally correlated) lognormal process ρ P (s) it is formally not lognormal. However, the sum of independent lognormal random variables is well approximated by a lognormal distribution with matched first and second moment (Fenton 1960;Marlow 1967). Here, we extend this approximation to the integral of temporally correlated lognormal processes. The first 2 moments of the neuron density follow from the averages of Equation (9), The lognormal approximation with matched moments is parameterized by where we used x = e μ+ 1 2 σ 2 and x 2 = x 2 [e σ 2 − 1] for a lognormal variable x. Note that the parameters of the lognormal distribution are the mean of the logarithmic neuron density μ N (t) = ln ρ N (t) and the variance of the logarithmic neuron density σ N (t) 2 = (ln ρ N (t)) 2 .

Across areas
Thus far, the model accounts for the lognormal distribution of neuron densities within an area. Across areas, we hypothesize that the distribution of proliferation times (Rakic 2002;Cadwell et al. 2019) is the most important cause of the variability in the neuron densities. To characterize an individual area, we consider the average density ρ N (t) . For convenience, we introduce the auxiliary quantity such that the logarithm of the mean neuron density simplifies to where we inserted Equation (7) into Equation (8) to obtain the mean of the progenitor cell density, which we then plugged into Equation (10). In order to arrive at a lognormal distribution of the mean density ρ N (t) across areas, the terms on the r.h.s. of Equation (13) have to be normally distributed. In particular, this means that the proliferation times have to be distributed such that y(t) is Gaussian. The distribution p(y) of y = y(t), which is a strictly monotonic transformation of the proliferation time t, is related to the distribution of proliferation times p(t) through (Van Kampen 2007) The first factor on the r.h.s. can be obtained from Equation (12) and the second factor is a Gaussian probability density with mean μ y and variance σ 2 y . Hence, Equation (14) fully specifies the distribution of proliferation times; conversely, for proliferation times distributed according to Equation (14), the mean neuron density ρ N (t) is lognormally distributed across areas. Note that since the Gaussian has support on the entire real line, the neuron density needs to diverge for t → ∞. Thus, we restrict λ P (t), λ N (t), and σ (t) such that the neuron density would diverge in the hypothetical limit of an infinite proliferation time.

Parameter estimation
Above, we showed how a noisy rate of progenitor proliferation leads to a lognormal distribution of progenitor cell densities and neuron densities within an area and, for the distribution of proliferation times (14), to a lognormal distribution of (within-area) mean neuron densities across areas. In order to compare the predictions of the model with the data, we estimate the model's parameters using the available experimental data.
We restrict the analysis to the simplified case of constant rate λ P , λ N , and noise intensity σ and identical rates of progenitor and neuron proliferation, λ P = λ N ≡ λ. In particular the assumption of a constant rate (note that this is not a necessary assumption for the above theory) is simplifying because the cell cycle length varies during development (Kornack and Rakic 1998). Despite this simplifying assumption, the model quantitatively matches the data for the parameters inferred below (Fig. 5).
First, we determine the rate λ from the average cell cycle length of progenitor cells determined from a short period of 2 hours (Kornack and Rakic 1998) such that f luctuations and the conversion of progenitor cells to neurons can be neglected. For a given cell cycle length , the number of cells increases as 2 t/ = exp t ln(2)/ . Thus, the cell cycle length corresponds to a proliferation rate λ = ln(2)/ . Using the average cell cycle length of ≈ 1.5 days from macaque (Kornack and Rakic 1998;Picco et al. 2018), we obtain λ ≈ 0.46 days −1 . The proliferation time of areas varies between 30 and 60 days in macaque (Rakic 2002) which we expect to be similar in the marmoset since macaques and marmosets have similar gestation times of 5.5 and 4.5 months, respectively (Schultz-Darken et al. 2016). Thus, we set the median proliferation time per area to t 1/2 = 45 days, which determines μ y = y(t 1/2 ) because the median of the distribution (14) is given N (y | μ y , σ 2 y )dy and the median of the normal distribution is at y = μ y .
In addition, the mean μ y of the auxiliary variable y is also constrained by the distribution of the variance σ 2 N of the logarithmic neuron density across areas; we will use this additional constraint to determine the noise intensity σ . For a fixed proliferation time t, σ N (t) 2 is given by Equation (11). Since σ 2 N is a strictly monotonically increasing function of t, its distribution can be written in terms of . Using Equation (14) for the distribution of proliferation times p(t) and Equation (12) to relate y(t) and t, the distribution of σ 2 N across areas is thus where f (y) = σ 2 N (t(y)) with t(y) the inverse of which follows from evaluating Equation (12) with λ P = λ N ≡ λ and constant σ 2 , and f −1 (σ 2 N ) is the inverse of f (y). Explicitly, it is given by which follows from inserting t(y) into σ 2 N (t) determined by Equation (11) and solving it for y. The maximum likelihood estimatorμ y for μ y with the likelihood p(σ 2 N ) is the empirical average of f −1 (σ 2 N ) across areas because the Jacobian 1/|f (f −1 (σ 2 N ))| in Equation (15) does not depend on the parameters and hence does not affect the Gaussian likelihood of μ y . Using the marmoset data, we obtain μ y ≈ 3.07. We finally arrive at σ ≈ 0.061 days −1/2 by enforcinĝ μ y = y(t 1/2 ), where y(t) is given by Equation (16), with t 1/2 = 45 days. With σ and μ y determined, we choose σ y such that less than 2% of the proliferation times are smaller than 30 days or larger than 60 days. This leads to σ 2 y ≈ 0.02. It remains to estimate the initial progenitor density ρ 0 . To this end, we consider the distribution of ln ρ N . For a fixed proliferation time t, ln ρ N (t) is determined by Equation (13). Thus, it is also a monotonic transformation of the random variable t, with distribution p(ln ρ N ) = dt d ln ρN p t(ln ρ N ) . Using the linear dependence between ln ρ N and y, Equation (13), in combination with the distribution of proliferation times (14), we obtain The maximum likelihood estimator for ln ρ 0 + μ y is the empirical average of ln ρ N across areas; subtracting μ y we obtain ρ 0 ≈ 3.8 × 10 3 cells/mm 3 . Note that in the simplified case considered here, the resulting distribution of proliferation times approaches a lognormal distribution in the left tail and a Gaussian in the right tail: for small times, 1 2 σ 2 t 1, a Taylor expansion of y(t) leads to y(t) ≈ ln(λt) and thus a lognormal distribution; for large times, 1 2 σ 2 t 1, y(t) grows linearly with y(t) ≈ ln 2λ σ 2 + 1 2 σ 2 t and thus the distribution approaches a Gaussian.
Simulation details We solve the SDE (5) using the Euler-Maruyama method with time step t = 0.05 in total N sample times with identical parameters for each of the 114 areas. Because the number of samples per area N sample varies, we randomly choose N sample for each area following N sample ∼ Poisson(36.6) where 36.6 is the average sample size per cortical area in the marmoset data. The subsequent analysis of the model data is identical to the analysis of the experimental data.

Lognormal distribution of neurons across cortical areas
We consider the neuron density distribution across cortex for several species (see Section Cell density data). The SW test (see Section Lognormality testing) concludes that the normality hypothesis of ln(ρ) cannot be rejected for mouse, marmoset, macaque 1 , human, galago 1 , owl monkey, and baboon (Fig. 1B). For the data sets macaque 2 and galago 2 , the normality hypothesis is rejected (P< 0.05); however, in these data sets, the densities were sampled neither uniformly nor based on a cytoarchitectonic parcellation. The normality hypothesis for the distribution of logarithmic densities across cytoarchitectonic areas is further supported by Figure 1C, which shows that the relation between theoretical quantiles and ordered samples is almost perfectly linear except for macaque 2 and galago 2 .
For lognormality testing, we removed the large outliers (marked with a red cross in Fig. 1C). The outliers are area V1 in macaque 1 and marmoset, which have densities far outside the range for all other areas in both species, and area APir in marmoset, which has a noticeably distinct cytoarchitecture with respect to the rest of the cerebral cortex (Atapour et al. 2019).
Next, we test the z-scored ln(ρ) from the different species and data sets against each other and find that they are pairwise statistically indistinguishable (α = 0.05 level; two-sample twosided Kolmogorov-Smirnov test, see Fig. S1 for full test results).
Additionally, we control for cell types in the distributions of the mouse, galago 1 , owl monkey, and baboon data. In the mouse data, different types of neurons and glia were labeled with specific genetic markers and their respective densities were reported separately for all cell types (Erö et al. 2018). In the galago 1 , owl monkey, and baboon data sets, the total numbers of cells and neurons were reported separately (Collins et al. 2010). We show that the neuron density distributions for all subtypes of neurons in the mouse are compatible with a lognormal distribution ( Fig. S2; SW test on ln(ρ), P > 0.05) while glia are not-with the notable exception of oligodendrocytes. When neurons and glia are pooled together ( Fig. S2C and D), the distribution of ln(ρ) still passes the SW normality test, likely due to the distribution being dominated by the neurons. Similar observations are made in the baboon data, where the glia do not pass the lognormality test, but the neurons do. In the cases of galago 1 and owl monkey both the neurons and glia pass the lognormality test (Fig. S2), which may, however, be partly due to the small number of density samples (N = 12 in both cases). Thus, the mouse and baboon data-with large samples sizes (N = 42 and N = 142, respectively)-suggest that it is the neuron densities that follow a lognormal distribution but not necessarily the glia densities. Finally, we perform a control test on the different staining types-Nissl and NeuN-using the macaque 1 data. The staining methods differ in their treatment of glia: NeuN tends to label neuronal cell bodies only while Nissl indiscriminately labels both neurons and glia (Yurt et al. 2018). We show that regardless of staining type the cell densities pass the lognormality test ( Fig. S3; SW test on ln(ρ) with P > 0.05), suggesting that counting some glia in the cell densities does not confound our analysis of the macaque 1 data.
Taken together, the normality test, the quantile plots, the pairwise tests, the cell-type comparison, and the staining method comparison provide compelling evidence that the logarithmized neuron densities are normally distributed across cytoarchitectonic areas. This also holds for uniformly sampled neuron densities (baboon) but not for a sampling that is neither uniform nor based on a cytoarchitectonic parcellation (macaque 2 , galago 2 ). Thus, the neuron densities are consistent with a lognormal distribution across the different cortical areas, as long as sampling is uniform.
The observation of a lognormal distribution across cortical areas raises the question whether there is a spatial pattern of densities across cortex consistent for all species. To address this question, we visualized the neuron or cell density over f lattened representations of the cortex for all data sets having a f lat map (8 out of the 9 included in this study; Fig. 2). Consistent with previous reports (Collins et al. 2010;Erö et al. 2018;Atapour et al. 2019;Beul and Hilgetag 2019), we observe a clear posterior-to-anterior density gradient in all the species shown in this study. Visual areas display the highest densities, whereas motor and frontal areas display the lowest density.

Lognormal distribution of neurons within marmoset cortical areas
To investigate whether the lognormal distribution holds within cortical areas, we leverage detailed estimates of neuron density in marmoset (Atapour et al. 2019). Neurons were counted within 150 × 150 μm counting frames for 4 strips per cortical area, all originating from the same subject. The within-area distributions of the sampled neuron densities ρ s across the counting frames in 3 representative areas (MIP, V2, and V3; Fig. 3A) again suggest a lognormal distribution. As before, we check for lognormality by testing ln(ρ s ) for normality with the SW-test (for full test results see Table S2). At significance level α = 0.05, the normality hypothesis is not rejected for 89 out of 116 areas; whereas at α = 0.001, this is the case for 114 out of 116 areas ( Fig. 3B and C). Thus, regardless of the precise significance threshold, the lognormality hypothesis cannot be rejected within most cortical areas in the marmoset cortex. The contribution to the neuron density distributions differs across cortical depth: the highest densities tend to occur in Layer 2 or around the center of the gray matter depth, whereas the lowest values appear either in the upper layers or near the white matter boundary (Fig. 3D). Taken together, the findings from Figs. 1 and 3 show that the lognormal distribution of neuron densities can be found at different scales, both within and across the cortical areas.

Statistical model comparison
To complement the statistical hypothesis tests on the logarithmic densities, we compared the lognormal model with 6 other statistical distributions based on the relative likelihood (see Section Statistical model comparison). We included statistical distributions with support in R + since neuron densities cannot be negative: lognormal, truncated normal, inverse normal, gamma, inverse gamma, Lévy, and Weibull. Of these, the lognormal, inverse normal, and inverse gamma distributions stand out as the distributions with the highest relative likelihoods, both across the entire cortex and within cortical areas (Fig. 4A, Fig. S4A). A visual inspection of the fitted distributions reveals that the lognormal, inverse normal, and inverse gamma distributions produce virtually indistinguishable probability densities (Fig. 4B,  Fig. S4C); indeed, the relative likelihoods of the 3 models are above 0.05 in all cases. This suggests that the data could theoretically be distributed according to either the lognormal, inverse normal, or inverse gamma distribution. To narrow down the model comparison, we show below how the lognormal distribution could arise from the simple biophysical process of noisy cell division. In contrast, we are not aware of a simple mechanism that could give rise to inverse normal or inverse gamma distributions.

A minimal model for the emergence of lognormally distributed neuron densities
The finding of lognormally distributed neuron densities raises the question how the intricate process of neurogenesis (Dehay and Kennedy 2007;Rakic 2009;Cadwell et al. 2019) culminates in this distribution within and across areas in several mammalian species.
On the within-area level, a minimal model shows that there is no need for a specific regulatory mechanism (see Model of neurogenesis with variability for further details): assuming that the proliferation of the neural progenitor cells is governed by a noisy rate λ P (t) + ξ(t), where λ P (t) denotes the mean rate and ξ(t) is a zero-mean Gaussian white noise, the resulting density of progenitor cells is lognormally distributed.
The cells produced by the cell division of the neural progenitor cells are terminally differentiated neurons, which thus do not divide further. This renders the final neuron density a linear integral of the (changing) neural progenitor cell density. While this additive process formally does not preserve the lognormality of the neural progenitor density, the resulting neuron density distribution is statistically indistinguishable from a lognormal distribution with matched moments (SW test P = 0.4 with N = 2, 000 samples in Fig. S5F)-akin to the lognormal approximation for the sum of independent lognormal random variables (Fenton 1960;Marlow 1967). Put differently, the lognormally distributed density of neural progenitor cells leads to an approximately lognormally distributed density of neurons. Thus, the lognormal neuron density distribution within areas could be a hallmark of a progenitor cell proliferation process with variability. On the across-area level, the proliferation times for each area become relevant, since they vary up to twofold (Rakic 2002). We therefore hypothesize that the distribution of proliferation times is the most important source of the variability across areas (variability due to area-specific proliferation rates (Lukaszewicz et al. 2005) is discussed below). Since the average neuron density within an area is determined by the proliferation time, the distribution of proliferation times specifies the distribution of area-averaged neuron densities across areas. This relation can be inverted to determine the distribution of proliferation times from the lognormal distribution of (within-area) average neuron densities (see Model of neurogenesis with variability). A specific prediction of this model is that the logarithmic mean ln( ρ ) and variance of ln(ρ) are related through the proliferation time-both mean and variability increase with proliferation time in approximate proportion to each other (see Equations (11) and (13), respectively). Indeed, we observe a linear correlation in the marmoset data (Pearson r = 0.51, P ≤ 10 −8 , Fig. 5C) as well as in the data produced by the model (Pearson r = 0.4, P ≤ 10 −4 , Fig. 5D).

Discussion
In this work, we show that neuron densities are compatible with a lognormal distribution across cortical areas in multiple mammalian cortices and within most cortical areas of the marmoset, uncovering a ubiquitous organizational principle of cerebral cortex. The distributions of neuron and cell densities in general depend on the underlying spatial sampling: We found that neuron densities follow a lognormal distribution within cytoarchitectonically defined areas, across such areas, and when averaged within small parcels uniformly sampled across cortex, but not when sampled in a highly non-uniform manner not following cytoarchitectonic boundaries.
Furthermore, we show that none of a sizeable list of statistical models outperform the lognormal distribution. Our results are in agreement with the observation that surprisingly many characteristics of the brain follow lognormal distributions (Buzsáki and Mizuseki 2014). Moreover, this analysis highlights the importance of characterizing the statistical distributions of brain data because simple summary statistics-such as the mean or standard deviation-lack nuance and are not necessarily a good representation of the underlying data.
These findings are based on 9 publicly available data sets for 7 different species. While the majority of these have sample sizes in the range 36-114, for completeness we also included 2 data sets  consisting of 10 samples each. As the latter sample size does not lead to a powerful test of lognormality, it would be desirable to repeat the analysis with more extensive data once available. In addition, data from multiple subjects of the same species would allow testing the consistency of the lognormal distribution across individuals.
Finally, we propose a minimal model that accounts for the emerging lognormal distributions based on a noisy cell division process of the neural progenitor cells and their specification to neurons. In principle, a multiplicative process with Gaussian noise leads to lognormal distributions at any scale; the additive specification from progenitors to neurons approximately preserves the lognormality. However, this noisy process alone cannot explain the coexistence of lognormal distributions of ρ at different scales (within and across areas). When the distributed proliferation times are considered, the model can account for both within-area and across-area distributions.
The model explains the observed lognormality based on a minimal set of assumptions; hence, we did not include further mechanisms like cell death (apoptosis), migration, volumetric growth, the generation of other postmitotic cells, or area-specific proliferation rates. However, none of these additional mechanisms affect the main conclusions. Apoptosis is a widespread phenomenon (Oppenheim 1991;Inglis-Broadgate et al. 2005;Kalinichenko and Matveeva 2008) during neurogenesis and it can be modeled as a multiplicative process affecting the neuron density alone. Thus, the neuron density still depends linearly on the progenitor cell density and approximately maintains lognormality. While apoptosis could add to the inter-area variability, our model shows that it is not necessary for the distribution of variability per se-in agreement with Oppenheim et al. (1989) who found that the final pattern of spatial variability in spinal motoneuron density was already present before the onset of cell death and migration. Following the radial unit hypothesis (Rakic 2009), we focus on radial migration of progenitor cells, ignoring tangential migration. In our model, this assumption can be relaxed as long as there is no net increase or decrease in progenitor cell density due to tangential migration. Furthermore, we modeled the cell density in a fixed final target volume; thus, the effects of volumetric growth do not affect the model. The generation of other postmitotic cell types would reduce the effective proliferation rate of progenitor cells; this only leads to quantitative but not qualitative changes in the resulting distributions. Finally, it has been shown that cortical areas can have different proliferation rates (Polleux et al. 1997;Lukaszewicz et al. 2005;Dehay and Kennedy 2007). If the proliferation rates are constant in time and lognormally distributed across areas, this additional variability would broaden the neuron density distribution but preserve the lognormal shape. However, in contrast to distributed proliferation times, it would not lead to the correlation between mean density and variability seen in the data (Fig. 5C), because the proliferation rate affects the mean density but not the variance of the logarithmic progenitor density (see Eq. (7)). Furthermore, to the best of our knowledge, a difference in proliferation rate has only been shown for areas V1 and V2 (Lukaszewicz et al. 2005; Dehay and Kennedy 2007)-thus, we speculate that this is an important reason for the drastically higher neuron density in V1.
In contrast to the neuron densities, we observed that the densities of most glia types are not compatible with a lognormal distribution (Fig. S2). Across brain areas, our model requires distributed proliferation times to explain the emergence of the lognormal distribution of neurons. However, it is established that gliogenesis occurs after neurogenesis in both rodents and humans (Miller and Gauthier 2007;Semple et al. 2013), and glia are likely to have proliferation times distinct from those of neurons. Thus, given different glia proliferation times, our model does not necessarily predict a lognormal distribution across cortical areas, in agreement with the statistical tests from experimental data (Fig. S2). The principles governing the within-area distributions should, however, be very similar for both glia and neurons; thus, within cortical areas, our model predicts lognormal distributions of glial densities. Unfortunately, we are not aware of any experimental quantitative data on which this prediction could be tested.
Given the fact that cortical space is limited, one may expect a negative correlation between neuron density and soma size. We investigated this for the human data set, the only data set providing soma sizes. Indeed, there is a significant correlation (Pearson r = 0.67, P ≤ 10 −30 , Fig. S6) between inverse soma sizes and cell density in the human data set. Here, soma sizes were measured in terms of mean surface areas from the Nissl-stained slices by approximating the soma shapes as ellipses (π/4 times height times width). Since the cell densities were obtained linearly from the 2D densities, the total surface area taken up by the cells is proportional to their density times their mean surface area. Despite the positive correlation, this product is far from constant (Fig. S6), and hence also the area taken up by the extracellular space, neurites, and other structures such as blood vessels varies across areas and layers. To disentangle how these various factors are related to neuron density, joint measurements of the different components would be needed.
In principle, cortex-wide organizational structures might be by-products of development or evolution that serve no computational function (Otopalik et al. 2017)-but the fact that we observe the same organizational principle for several species and across most cortical areas suggests that the lognormal distribution serves some purpose. Heterogeneous neuron densities could assist computation through their association with heterogeneity in other properties such as connectivity and neuronal time constants (Rall 1969;Hilgetag et al. 2019); indeed, such heterogeneity is known to be a valuable asset for neural computation (Duarte and Morrison 2019;Perez-Nieves et al. 2021). Alternatively, localized concentration of neurons in certain areas and regions could also serve a metabolic purpose since centralization might support efficient transport of metabolites across neurons and astrocytes (Bélanger et al. 2011;Magistretti and Allaman 2015). Energy efficiency is particularly relevant since a large portion of the brain's energy consumption is used to support the communication between neurons (Attwell and Laughlin 2001;Laughlin and Sejnowski 2003). Also from the perspective of cortical hierarchies it makes sense to have few areas with high neuron densities and many areas with lower neuron densities: Low-density areas contain neurons with large dendritic trees (Elston and Rosa 1998) receiving convergent inputs from many neurons in high-density areas lower in the hierarchy. The neurons with extensive dendritic trees in higher areas are involved in different, area-specific abstractions of the low-level sensory information (Kumar et al. 2007;Brincat et al. 2018). All in all, there is probably not a single factor that leads to lognormal neuron densities in the cortex; further research will be needed to refine our findings and uncover the functional implications. Conf lict of interest statement. None declared.

Data and materials availability
This work produced no new data and instead relied on a corpus of neuron density data available from the literature, which we gratefully acknowledge; see Section Cell density data for a detailed description. The code used to analyze the data and create all figures is available in an open repository (DOI: 10.5281/zenodo.7871536).